estimate<-F
use_GA<-F
use_kappa_instead_of_theta<-F
num_max_iterations<-1
verbosity<-F

   source("model.R")

## for estimation, set A_array to the one read in from data
  A_array_months<-A_array_months_data

calculate_log_likelihood<-function(estimation_parameters_in_vec)
	{
  # set parameters to those read in from file
  parameters_in_vec<-parameters_readin$estimate 
  
  #update with those in optimizer
  parameters_in_vec[!parameters_readin$fixed]<-estimation_parameters_in_vec 
  if (verbosity==T) {
    print(parameters_in_vec[!parameters_readin$fixed])
  }
  
  parameters_in<-relist(parameters_in_vec, parameters)
  
  n_parameters_out_of_bounds<-sum(rowSums(check_pars(parameters_readin)(parameters_in)))
  
  if (n_parameters_out_of_bounds>0)  {
    list(neg_log_like = n_parameters_out_of_bounds*20000)
  } else {
    
  attach(parameters_in)
  
	mu_vec_y<-with(student_data, hsgpa*omega_vec[1] + combact*omega_vec[2] + black*omega_vec[3] + male*omega_vec[4] + studyhs*omega_vec[5] + estudy*omega_vec[6])
  mu_vec_study<-pmax(with(student_data, hsgpa*omega_vec[7] + combact*omega_vec[8] + black*omega_vec[9] + male*omega_vec[10]                           + studyhs*omega_vec[11] + estudy*omega_vec[12]), -10)
  
  someone_has_zero_study<-F
  
	### Solve model for each month
	for (month in 1:n_months)
		{
		soln<-calc_equilibrium(A_array_months[,,month])(mu_vec_study) 

		s_array_months[,month]<-soln$ans

		s_not_i_array_months[,month]<-calc_s_not_own(A_array_months[,,month])(s_array_months[,month])
		}

  if(min(s_array_months, na.rm=T)==0)
  {
    detach(parameters_in)
    print("someone has zero studytime")
    list(neg_log_like = 20000)
    
  } else
    {
      
  ## takes model study times and maps to studytime likelihood
  	for (survey in 1:n_study_surveys)
    	{
  		likelihood_s_array[,survey]<-Vectorize(likelihood_s)(s_obs=student_data[,paste0("s", survey)], s_model=s_array_months[,month_corresponding_to_survey[survey]])
      ## ignore students with 0 friends
        if (ignore_0_friend_students==T){
          likelihood_s_array[rowSums(A_array_months[,,month_corresponding_to_survey[survey]])==0,survey]<-1     
        }
      }
  
    for (semester in 1:2)
  		{
  		## average study time per semester
        s_i_avg_array[,semester]<-s_array_months[,semester]
      
  		## average friend study time per semester
        s_not_i_avg_array[,semester]<-s_not_i_array_months[,semester]
      
  		  y_array[,semester]<-y(s_i_avg=s_i_avg_array[,semester],s_not_i_avg=s_not_i_avg_array[,semester], mu_i=mu_vec_y)
  		  likelihood_y_array[,semester]<-Vectorize(likelihood_y)(y_obs=student_data[,paste0("gpa", semester)], y_model=y_array[,semester])
      
  		  y_growth_array[,semester]<-y(s_i_avg=s_i_avg_array[,semester],s_not_i_avg=s_not_i_avg_array[,semester], mu_i=0)
      
      ## ignore students with 0 friends?
  		if (ignore_0_friend_students==T){
    		  likelihood_y_array[rowSums(A_array_months[,,semester])==0,semester]<-1     
    		}
      }
  
  	sum_log_likelihood<-sum(log(likelihood_s_array), na.rm=T) + sum(log(likelihood_y_array), na.rm=T)
  
  	detach(parameters_in)
  
  	list(neg_log_like = -sum_log_likelihood, pars = parameters_in, s_array_months_model = s_array_months, y_model_array = y_array, mu_y_model=mu_vec_y, mu_study_model = mu_vec_study, likelihood_s_array=likelihood_s_array, likelihood_y_array=likelihood_y_array, y_growth_array=y_growth_array)
    }
  }
}

for_optim<-function(parameters_in) calculate_log_likelihood(parameters_in)$neg_log_like
for_GA<-function(parameters_in) -calculate_log_likelihood(parameters_in)$neg_log_like

## helper functions
  update_theta_given_kappa<-function(par_df)
  {
    out_df<-par_df
    
    kappa_vec_tmp<-with(out_df, out_df[par_name %in% paste0("kappa_vec",1:4),])$estimate
    theta_vec_tmp<-vector(length=4)
    theta_vec_tmp[1]<-with(out_df, out_df[par_name %in% "beta_vec2",]$estimate) - kappa_vec_tmp[3]
    theta_vec_tmp[2]<- -kappa_vec_tmp[4]
    theta_vec_tmp[3]<- -kappa_vec_tmp[1]
    theta_vec_tmp[4]<- -kappa_vec_tmp[2]
    out_df[out_df$par_name %in% paste0("theta_vec",1:4),]$estimate<- theta_vec_tmp
    
    out_df
  }  
  
  update_kappa_given_theta<-function(par_df)
  {
    out_df<-par_df
    
    theta_vec_tmp<-with(out_df, out_df[par_name %in% paste0("theta_vec",1:4),])$estimate
    
    kappa_vec_tmp<-vector(length=4)
    kappa_vec_tmp[1]<- -theta_vec_tmp[3]              
    kappa_vec_tmp[2]<- -theta_vec_tmp[4]
    kappa_vec_tmp[3]<-  with(out_df, out_df[par_name %in% "beta_vec2",]$estimate) - theta_vec_tmp[1]
    kappa_vec_tmp[4]<- -theta_vec_tmp[2]
    
    out_df[out_df$par_name %in% paste0("kappa_vec",1:4),]$estimate<- kappa_vec_tmp
    
    out_df
  }
  
  compare_pars<-function(oldpars, newpars){
    changed_pars<-newpars$fixed==F
    data.frame(par_name=oldpars$par_name,old=oldpars$estimate,new=newpars$estimate)[changed_pars,]
  }

estimating_these_pars<-parameters_readin[parameters_readin$fixed==F,]

if(estimate==TRUE){
  parameters_old<-parameters_readin
  guess_vec<-parameters_readin$estimate[!parameters_readin$fixed]
  
	# don't use genetic algorithm
  if(use_GA==F){
    num_max<-0 # number of maximization iterations
    
    while (num_max<num_max_iterations){
      
      ans<-optim(guess_vec, for_optim, method="Nelder-Mead", control=list(trace=6, maxit=15000, parscale=parameters_readin$scale[!parameters_readin$fixed], alpha=1, beta=0.5, gamma=1.5))
      save(ans, file="./estimation_stuff/ans")
      parameters_writeout<-parameters_readin
      parameters_writeout$estimate[!parameters_readin$fixed]<-ans$par

      if(use_kappa_instead_of_theta==F)
      { # using thetas to estimate, so update kappa
        parameters_writeout<-update_kappa_given_theta(parameters_writeout)
      } else {
        # using kappas to estimate, so update theta
        parameters_writeout<-update_theta_given_kappa(parameters_writeout)
      }
    
      write.csv(parameters_writeout, file="./estimation_stuff/parameters_writeout.csv", row.names=F)
      parameters_readin<-parameters_writeout
      guess_vec<-parameters_readin$estimate[!parameters_readin$fixed]
      num_max<-num_max+1
      print(num_max)
      print(compare_pars(parameters_old, parameters_readin))
    }
  }
  
  # use genetic algorithm
  if(use_GA==T){
    GA_min<-pmax(parameters_readin$lower_bound+1/100, -1)
    GA_max<-pmin(parameters_readin$upper_bound-1/100, 3)
    GA_min<-pmax(GA_min, parameters_readin$estimate - abs(parameters_readin$estimate)*0.5)
    GA_max<-pmin(GA_max, parameters_readin$estimate + abs(parameters_readin$estimate)*1.5)
    GA_min<-GA_min[!parameters_readin$fixed]
    GA_max<-GA_max[!parameters_readin$fixed]
    
    GA<-ga(type="real-valued", fitness=for_GA, min=GA_min, max=GA_max, popSize=20, pmutation=0.1, maxiter=50, suggestions=guess_vec)
		GA_ans<-summary(GA)
		save(GA_ans, file="./estimation_stuff/GA_ans")
	  parameters_writeout<-parameters_readin
	  parameters_writeout$estimate[!parameters_readin$fixed]<-GA_ans$solution
    
    if(use_kappa_instead_of_theta==F)
    { # using thetas to estimate, so update kappa
      parameters_writeout<-update_kappa_given_theta(parameters_writeout)
    } else {
      # using kappas to estimate, so update theta
      parameters_writeout<-update_theta_given_kappa(parameters_writeout)
    }
    
	  write.csv(parameters_writeout, file="./estimation_stuff/parameters_writeout.csv", row.names=F)
	  parameters_readin<-parameters_writeout

	  print(compare_pars(parameters_old, parameters_readin))
  }
} else {
  print(paste("likelihood at parameters_readin:", for_optim(parameters_readin$estimate[!parameters_readin$fixed])))
}

  parameters_new<-relist(parameters_writeout$estimate, parameters)

  
